Effects of spin vacancies on magnetic properties of the Kitaev-Heisenberg model 
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We study the ground state properties of the Kitaev-Heisenberg model in a magnetic field and 
explore the evolution of spin correlations in the presence of non-magnetic vacancies. By means 
of exact diagonalizations, the phase diagram without vacancies is determined as a function of the 
magnetic field and the ratio between Kitaev and Heisenberg interactions. We show that in the 
(antiferromagnetic) stripe ordered phase the static susceptibility and its anisotropy can be described 
by a spin canting mechanism. This accounts as well for the transition to the polarized phase when 
including quantum fluctuations perturbatively. Effects of spin vacancies depend sensitively on the 
type of the ground state. In the liquid phase, the magnetization pattern around a single vacancy in 
a small field is determined, and its spatial anisotropy is related to that of non-zero further neighbor 
correlations induced by the field and/or Heisenberg interactions. In the stripe phase, the joint effect 
of a vacancy and a small field breaks the six-fold symmetry of the model and stabilizes a particular 
stripe pattern. Similar symmetry-breaking effects occur even at zero field due to effective interactions 
between vacancies. This selection mechanism and intrinsic randomness of vacancy positions may 
lead to spin-glass behavior. 

PACS numbers: 75.10.Jm, 75.30.Et, 71.55.-i, 71.30.-|-h 
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I. INTRODUCTION 

Much attention has been paid recently to the existence 
and properties of spin liquids in strongly correlated ma- 
terials, both experimentally and theoreticalljiii^. Such 
phases, first proposed by Anderson, have been put for- 
ward in the context of high Tc superconductivity as a 
possible precursor state, but they can occur in a much 
broader class of situations, both in real materials^"— and 
in microscopic models2r2£. Their characterization is chal- 
lenging - on the experimental side, it requires generally 
involved low-temperature techniques to prove the ab- 
sence of long-range magnetic order down to temperature 
T = 0, and concerning the theory, either in Hubbard-likc 
models or in localized spin models, spin liquids usually 
result from frustrating interactions but compete with a 
variety of possible ordered phases, from usual antiferro- 
magnetic (AF) order to valence bond crystals, spin ne- 
matics etc.— 

Clear evidence of a spin liquid ground state was found 
from an exact solution of the Kitaev modelii of spins 
1/2 on the honeycomb lattice. This liquid ground state 
is characterized by a gapless spectrum (with low-energy 
excitations described in terms of Majorana fcrmions) and 
gapped vortices. The Kitaev model itself triggered a 
large variety of studies, focusing on topological prop- 
erties, (non)-abelian excitations, extensions to higher 
dimensionsi^, effects of magnetioi^ or non-magnetic im- 
purities and of a magnetic fiel d^^i^^" — . 

Besides its mathematical beauty and the original prop- 
erties mentioned above, the Kitaev model was recently 
foundii to be relevant to orbitally degenerate systems 
with strong spin-orbit coupling such as layered iridates 
Na2lr03 and Li2lr03. There, around each Ir atom, the 
octahedral environment [see Fig. HJa)] of oxygens results 
in a configuration with one hole in the t2g manifold; 



and the large spin-orbit coupling selects locally a low- 
energy doublet of states defining a pseudospin 1/2; and 
finally, superexchange processes between Ir ions driven 
by moderately strong Coulomb interactions lead to an ef- 
fective description by the Kitaev-Heisenberg mode l^^i^^ 
for pscudospins 1/2 residing on the honeycomb lattice. 
The study of this model showed the robustness of the 
spin liquid phase for finite AF Heisenberg interactions 
[small compared to ferromagnetic (FM) Kitaev interac- 
tions] , and a quantum phase transition between this liq- 
uid phase and a stripe-ordered phase. In Na2lr03, where 
interaction parameters could be close to this transition, 
AF order was found by magnetization and specific heat 
mcasurement o^^i^° and x-ray magnetic scattering^!. 

The question of how vacancies influence the order and 
the magnetic response in the Kitaev-Heisenberg model 
is of both experimental and theoretical interest. Indeed 
the vacancy can act as a probe allowing to measure spin 
correlations in its vicinity; and it can also have dras- 
tic effects on these correlations. For concrctcncss, in 
SU{2) spin-1/2 antifcrromagncts with long-range order, 
a non-magnetic impurity can enhance the staggered mag- 
netization around the vacanc y^^i^^ . In dimerized phases, 
the impurity is accompanied by a spinon, which is typi- 
cally confined to the impurity in Valence Bond Crystals 
and deconfined in Resonating Valence Bond liquidsSi, al- 
though the question of the confinement length can be 
subtle^^. In a Neel-ordered phase, the impurity results in 
a Curie- type contribution to the finite-T susceptibilitj*^^. 
In the spin-liquid phase of the Kitaev model, a magnetic 
moment is induced around a spinless vacancy; interest- 
ingly, its contribution to the spin susceptibility diverges 
logarithmically for T — )- 0^. This makes it relevant to 
investigate the vacancy-induced magnetic response in a 
model interpolating between SU{2) and Kitaev-like inter- 
actions. Moreover, in the iridate samples studied up to 
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now, a substantial amount of disorder might be present, 
e.g., site- mixing effects witli non- magnetic alkaline ions 
on the hexagonal Ir-sublattice^S, which should be taken 
into account for proper understanding of their physical 
properties. 

Our aim is to study the effect of non-magnetic vacan- 
cies (and of pairs of vacancies) in the different magnetic 
phases of the Kitaev-Heisenberg model, at zero temper- 
ature. We adopt similar notations as in Ref. [l^ with 
a parameter < a < 1 interpolating between Heisen- 
berg and Kitaev models, and a uniform magnetic field 

h — {hx, hy,hz) with amplitude h = yjh'^ + + ^z- 
H = -2a ^ a.>J + (l-«)^a,.a,-^/J.a„ (1) 

where 7 = x,y, z labels simultaneously an axis in spin 
space, and (in the interaction terms) a bond direction 
of the honeycomb lattice. For convenience we express 
H in terms of Pauli matrices ct^. We take as unit cell 
a z-hond [see Fig. [llb,c)], on which Kitaev interactions 
are of the form crlo'j- Each site is labeled by an index 
«, in bijection with (r, /3) where r denotes the position of 
the unit cell and /3 = A, B the sublattice index, respec- 
tively. Elementary translations between neighbor cells 
are ni/2 = (±1/2, •\/3/2) in Cartesian coordinates. We 
consider periodic clusters of = 24 sites (this cluster 
can be represented with an hexagonal shape, see Fig.[n|), 
iV = 32 or = 16 sites (both with the shape of a paral- 
lelogram). These clusters respect all spatial symmetries 
of the model, i.e. translation and 7r-rotation - other sym- 
metries combining spin- and spatial rotations^^ are not 
used here. 

In Sectionini we discuss several features of the model in 
a magnetic field. We will show how a spin canting mech- 
anism allows one to understand the phase diagram for a 
field oriented along an easy axis, as well as to understand 
a directional anisotropy of the static susceptibility. Also 
the perturbations due to Hcisenberg interactions and due 
to small magnetic fields on the spin correlation functions 
are explored. In Section lllll we study the effect of a single 
non-magnetic vacancy in the system, in various phases of 
the model, and analyze the magnetization pattern in the 
vicinity of this vacancy in a small field. In Section IIVI 
we discuss situations where two vacancies, close enough 
to each other, effectively interact. We focus mainly on 
the stripe phase and explain the underlying mechanism 
of local selection of an ordered pattern by a vacancy pair. 
Section Ivl provides a short summary and some conclud- 
ing discussions. 

II. VACANCY-FREE SYSTEM: SPIN 
CORRELATIONS, MAGNETIZATION, PHASE 
DIAGRAM 

We first consider some basic properties of the Kitaev- 
Heisenberg model in a magnetic field, and study the mag- 



netic response to a field oriented along different crystal- 
lographic directions. This comparison is important be- 
cause of the symmetry of the model, which is clearly not 
SU (2) invariant due to the Kitaev interactions. In the 
context of iridates, this anisotropy is due to spin-orbit 
coupling which results in the easy spin axes x || [100], 
y II [010], z II [001] [corresponding to the octahedral axes, 
see Fig. [TJa)], compared to the cubic diagonal [111] di- 
rection. The clearest anisotropic features are seen in the 
stripe-ordered phase, and we compare this situation to 
that of the Neel phase. 

A. The competing phases in a magnetic field 

We first want to characterize the various phases found 
in a magnetic field. It is already knowni^ that at zero 
field, the spin liquid phase extends from otsjj^ ~ 0.80 
to a = 1 (the Kitaev limit); the stripe-ordered phase, 
with spins pointing either in x, y or z direction [for the 
last case, the spatial structure, i.e. z-stripe pattern, is 
shown in Fig. ^^a)], is located between a^r/s — 0.40 and 
OLsjL\ and a Neel-ordered phase is found for < a < 
otNjs- To determine whether these phases persist in a 
magnetic field and to understand the magnetic response 
of the system, one needs adapted structure factors which 
act as order parameters. For the Neel phase, it reads: 

SNeel = 'Yp.''^^ei,li''\^r,e, ' ^rfi') (2) 
r.r' 

with v^^fii being 1 for /3 = /3', and —1 otherwise; for the 
stripe-ordered phase, the structure factors corresponding 
to X- and z- stripe patterns are defined as follows: 

with Qx/y/z shown in Fig. [Ijd); for instance, Qz = 
(0, ■^) which is equivalent to the more intuitive Q'^ ~ 
(27r,0) (but which is not in the first Brillouin zone). 
Both wave vectors account for a phase factor —1 be- 
tween neighboring unit cells separated by (±1/2, ±-\/3/2) 
in real space. These structure factors also evidence the 
symmetries spontaneously broken in each of the 6-fold 
degenerate ground states in the stripe phase: indeed at 
zero field, in addition to translational symmetries, the 
Hamiltonian possesses a Z2 * Z3 symmetry - Z2 stand- 
ing for the (Ji — > —Ui transformation and Z3 (labeled as 
C3 in Ref. [28l ) for the invariancc by a cyclic permutation 
of the X, y and z spin components coupled to a spatial 
rotation by 27r/3. 

Both AF (Neel and stripe) ordered phases of the model 
persist under moderate fields, which tilt spins by an an- 
gle proportional to the small field amplitude ( canted Neel 



3 










z 


B 





FIG. 1: (color online) (a) Schematic picture of an octahedron 
(6 oxygens surrounding a Ir*^ ion) with the octahedral axes x, 
y, z - which are easy axes on identically-labeled bonds in the 
Kitaev-Heisenberg model, (b) A hexagon in the Ir*""" plane, 
indicating bonds labeled x, y, z according to their local easy 
axis regarding Kitaev interactions, (c) Honeycomb (direct) 
lattice with the labeling x,y,z of bonds according to their 
in-plane orientation. A unit cell (centered on a z-bond) is 
represented in the dashed rectangle, as well as elementary 
translations ni/2. (d) First Brillouin zone, with the position 
of stripe-order wave vectors (7 = x, y, z). 
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FIG. 2: (color online) (a) Representation of a jz-stripe pattern 
(stripe pattern with spins oriented along the z = [001] axis 
- neighboring spins are aligned on z bonds and antialigned 
on other bonds, (b) Schematic representation, at the level of 
a single hexagon, of a canted stripe pattern: here a z-stripe 
pattern, with spins tilted of an angle 9 towards a field along 
[100] direction. 



the thermodynamic limit leads to an anisotropic zero- 
field susceptibility. On the other hand, the ground state 
of a finite cluster is known to be a singlet, which does not 
break SU{2) symmetry. In order to perform a classical 
treatment of the Heisenberg model in absence of sponta- 
neous symmetry breaking, one can consider for instance 
a SU{2) symmetric combination of classical Necl states - 
such a state can mimic, from an experimental viewpoint, 
the response of a polycristalline sample. The magnetic 
response of this state is isotropic in the sense that a small 
field results in a magnetization of amplitude mo along 
the field, and mg does not depend on the field direction. 
When considering the Kitaev-Heisenberg model for small 
finite a, with the classical approach above, one can find 
for any field direction a Neel pattern with spins perpen- 
dicular to it; within this approximation, the magnetic 
response is again isotropic, despite the absence of SU{2) 
symmetry due to small Kitaev interactions. However, 
the combination of quantum fluctuations and small Ki- 
taev interactions have been found to select one of the 
cubic axes as an easy spin axis and open thereby a small 
spin gapi^; this will lead to small anisotropies in the spin 
susceptibility. 

In the stripe phase, even at the classical level only six 
ordered patterns are allowed at zero field. For a generic 
field h with h^, hy and hz taking three distinct values, the 
three allowed stripe orientations are energetically differ- 
entiated, depending on the field direction. For instance, 
in a small field h \\ [100], the two z-stripe patterns are 
favored (along with the y-stripc patterns), since they 
allow for the field to tilt spins by a small angle 6 [see 
Fig.[5Jb)], such that the uniform magnetization develops 
an x-component. This effect is evidenced by structure 
factors characteristic of the ordered phase - in the case 
of Fig. [3] (with a = 0.6), in a small [001]-oriented field 
the x-stripc pattern becomes more favorable than the z- 
stripe pattern: S'^{Qx) jumps from approximatively 1/3 
at /i = to 1/2 at small finite h^, while S^{Qz) (not 
shown) drops to almost 0. As for the magnetization per 
site mo{hz), it grows linearly with hz in this regime. 

The tilting angle 9 can be determined simply in the 
classical approach, where the energy per unit cell as func- 
tion of 9, for given a and hz, is: 

Es{hz,e) = -(l + a) -2hzsme + 4:{l-a)sm'^e, (5) 



phased or canted stripe phased). To allow for this tilt- 
ing, the untiltcd spins (for h infinitesimal) should not be 
coUinear with the field^S, and the staggered magnetiza- 
tion of the selected ordered pattern has an angle with h 
as close to n/2 as possible. 

We first discuss the effect of a small field in the Neel 
phase, which allows to precise what we will understand in 
the following by (an)isotropy in the magnetic response. 
In the pure Heisenberg model, both classical and quan- 
tum cases, spontaneous breaking of SU (2) symmetry in 



and its minimization w.r.t. 9 gives a magnetization per 
spin Too = sin 00 = ^ila) ' relatively good agreement 
with the linear behavior observed at small field (see Fig. [3] 
for the case a — 0.6); also, the predicted a-dependence 
of the (uniform) susceptibility per spin, 

^ \dhz Jh=o 4(1 -a)' ^ ' 

is in good agreement with the numerical results in Fig. [S] 
Note that the direction of the field is relevant in the stripe 
phase: a field h \\ [001] differentiates stripe patterns. 
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FIG. 3: (color online) Magnetization mo (per site) and 
x— stripe structure factor S^{Qx) versus for various values 
of a and A'^. Dotted lines: Perturbative results for mo in the 
polarized phase. Inset: Phase diagram in the {a,hz) plane; 
continuous transition lines are interpolations of points from 
N = 24 data and dashed lines are perturbative estimates. 





' 1 


1 1 1 


- 


x-^ X (N=32) 






ao4(l-a)x' (N=32) 






(N=24) 




_ 


D-Q^fl'"! (N=32) 




o- 


e-- 




















V 

7 N. 






\ 

\ 






\ 

N 






S 

\ 

\ 


5 


0.6 0.7 


0.8 0.9 



a 

FIG. 5: (color online) Susceptibility per spin evaluated 
as mo/Ziz for small — 0.01, and its equivalent x'^^^' for a 
[lll]-oriented field, on clusters = 24 or iV = 32. Shown 
also the quantity 4(1 — a)x^ for = 32. 




FIG. 4: (color online) Evolution of correlations \{ala~j)\ with 
a on a A*' = 24 cluster at = 0, for various relative positions 
(i—j) with i and j indicated in the inset. The red dot indicates 
the value ~ 0.53 of the (0 — 1) correlation at a = 1, which 
coincides with an analytical result from Ref. [T^ . 



and selects two of them where spins are exactly orthog- 
onal to the field as /i — !• ; in contrast, in a /i |j [111] 
field the three stripe directions remain equivalent (thus 
S'^(Q^) = S'^(Q^) ~ 1/3 for /i 0); also, none of the 
uncanted stripe patterns is orthogonal to a [lll]-field, so 
the magnetic response for this field direction necessarily 
differs from the response to a [001]-field. 

We now turn to the spin liquid regime, corresponding 
to the weakly perturbed Kitaev model, and describe the 
effect of perturbations to the system. Our main focus 
is on the (cf (t|) correlations. In the pure Kitaev model 
(a = 1, /i = 0), these correlations are non-zero only 



when i,j are nearest neighbors on a z-bondii*i^, where 
they take a value (erf (T|)i"~^'''~"'' ~ 0.53 which wc can 
see in Fig. |4] (the finite-size correction to this quantity, 
i.e. to the ground state energy, is negligible even for the 
= 24 sites cluster). 

The gapless nature of excitations persists in a small 
axial field hz; in contrast, a [lll]-oriented field opens a 
gap resulting in a topologically ordered phaseiii^. Yet 
even in the first case, the field causes peculiar {afaj) 
correlations with spatial oscillations perpendicular to z- 
bonds and a power- law decaying envelope, as shown for 
a = 1 in Ref. Il6l Some of these trends appear in the spin- 
spin correlations displayed in Fig.[6l^a). Due to the small 
size, one cannot identify a power-law behavior; but one 
notices that correlations decrease with distance, much 
more slowly in the direction parallel to z bonds than 
in the direction perpendicular to them (i.e. parallel to 
fii — n2)- Given the periodic boundary conditions on 
the A^ = 24 lattice, this is in agreement with the pe- 
riod of three unit cells for spatial oscillations predicted in 
Ref. [l^. The liquid phase is stable with respect to small 
Heisenberg interactions (for a > ag/L — 0.80 at zero 
fieldi^). These interactions induce correlations between 
further neighbors, of different type than those induced 
by a small field [compare Fig.[SJa) and (b)]. The spatial 
alternation of sign on Fig. IHl^b) reminds of the z-stripe 
pattern, with the difference that here these correlations 
are short-ran ged . A proof of this last point has been 
given in Ref. [2^ for ft, = 0, and indeed absolute values 
[(cr^CT^)! of correlations shown in Fig. |4] seem to decrease 
rapidly - possibly exponentially ~ with distance, at least 
in the spatial direction orthogonal to z bonds. On near- 
est neighbors, these correlations are proportional to 1 — a 
(not too close to ag/j^). With both perturbations at play, 
i.e. when both 1 — a and arc finite and small, their 
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effects should add up; therefore, and since the correla- 
tions induced by Heisenberg terms decay exponentially, 
one should observe a correlation pattern with a power- 
law tail as in the case of a = ^ 0, but modulated 
as in a z-stripe pattern for short distances, at least for 
/i^ < 1 - a. 

Eventually, for any value of a, a sufficiently high field 
drives the system into a polarized phase with spins par- 
allel to h. In the case of a [001]-field, by considering 
transverse interactions erf crif and crfcr^ as perturbations 
[the unperturbed ground state being fully polarized for 
hz > 2(1 — a)], within second order perturbation theory 
one finds 



Too 



1 - 



4a - 2) 



(7) 



in good agreement with numerical data for > 1 (see 
dashed curves in Fig. [3]). 



B. Transitions between liquid, stripe-ordered, and 
polarized phases 

When a magnetic field hz is applied to a system which 
is initially stripe- or Neel ordered, aside from the spin 
canting in the AF ordered state, the excitation energy 
of a competing polarized state decreases with so that 
this state becomes the ground state above a critical field 
hg/p. This field can be evaluated by estimating the en- 
ergies (both per unit cell) Es{hz) of the canted stripe 
phase and Epzihz) of the polarized phase, as a function 
of fi.eld - then the equation Eg — Epz determines hg/p. 
For both quantities, we add to the classical value a term 
(last terms in both lines of Eq.[8]) accounting for quantum 
fluctuations, considered in second order perturbation the- 
ory, and evaluated respectively from the uncanted stripe- 
or the fully polarized state: 



(l-2a)2 



Epz{hz) = -2hz + 3-5a 



4 1 — a a 
2a^ 



h^+Aa-2 



(8) 



We compare this estimate of hg/p with numerical es- 
timates obtained from the /iz-dependence of the stripe 
structure factor S^{Qx) and the magnetization toq for 
iV = 24 and N = 32 clusters - both quantities show very 
weak size-dependence (see Fig. ^ and their derivatives 
with respect to hz are maximized for values very close 
to each other. For a = 0.6, this numerical estimation 
gives hg/p = 0.65(5) compared to a value 0.55 found 
from the criterion Eg — Epz- Both types of estimates 
are shown as a function of a, in the phase diagram in 
the inset of Fig. [31 For a > otg/h^ the magnetization 
curve niQ^hz) on a finite cluster at fixed a {a = 0.9 in 
Fig. [3]) shows a sharp increase in a field range centered on 
Hl/p ~ 0.15(1); this value is used to locate the transition 
to the polarized phase (our precision on hj^/p is limited 



here by larger finite size effects in the liquid state than 
in ordered phases). This can be understood by modeling 
the energy per unit cell for the liquid phase by 



Er 



3(1 - 3a)(f7 



n 



X 



(9) 



and looking for the field value such that E^ = Epz- This 
gives a value h^p ~ 0.13(3), which increases slightly with 
\ — a. Similarly, the position of the stripe-liquid transi- 
tion, which at zero field occurs at ag/i^, is almost insen- 
sitive to a small field hz- 

The schematic phase diagram in Fig. [3] has a similar 
structure as its counterpart for a [lll]-oriented field^^, 
except that we do not see a multicritical point at a = 
ag/L,hz ~ 0. This may be connected with the absence 
of C| symmetry in the Hamiltonian Eq.[l]in a [001]-field; 
moreover, the nature of the spin liquid for a > ag/j^ in a 
small axial field differs from that of the topologically or- 
dered phase in a small [lll]-field, mainly in the structure 
of the low-energy spectrum. 



C. Susceptibility anisotropy 

We analyze now in more detail the zero-field suscep- 
tibility: it has unusual properties due to a special spin 
symmetry of the model, which is lower than SU{2). In 
the stripe phase (ajv/s < o: < ag/Lj^ if we consider the 
general case of a [mnp]-oriented field, that is h parallel to 
rax + ny+pz, the corresponding susceptibility be- 



haves as 



, where Cm„p is a rational number which can 



be found from the classical picture of a spin-flop mech- 
anism. One can see in Fig. [S] that susceptibilities along 
different axes follow well this type of law in the stripe 
phase. The coefficient Cmnp is maximized for a [001]-field 
(cooi = 1/4, see Eq. |6]), but also for a [101]-field, and 
minimized for a [lll]-field (cm = 1/6). Note here that 
the [111] axis corresponds, in the context of susceptibility 
measurements in Ref. [13, to the c axis of Na2lr03, i.e. 
perpendicular to the Ir-honeycomb layers [see Fig.[ljb)]. 
Concerning the susceptibility in the ab plane, this model 
predicts that - in the stripe phase - x takes a value 



between 



Trr—, and x'"'' 

4(1 — q) 



An av- 



5(l-a)- 

erage over the field orientations in the ab plane gives 
4(l-a)(x''') = 2/3+-^ -0.94. 

No such simple expressions for the susceptibility are 
available for the liquid phase; but there, despite large 
finite-size effects (compared to the stripe phase) one no- 
tices that for a given axis susceptibilities tend to increase 
with a, with again x'^^^' 1^ x^™^'- ^s for the Neel phase, 
a classical approach does not predict there any anisotropy 
in the static susceptibility (see related discussion in Sec- 
tion |lT^, and the anisotropy resulting from quantum 
fluctuations at small finite a should be smaller than in 
the stripe phase - but the Neel phase is obtained for 
values of a likely unrealistic for Na2lr03i^ and also not 
supported by experiment^. 
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The experimental data of Ref. [g^ reveal an anisotropy 
in the static susceptibility as well, but the measured sus- 
ceptibility is significantly larger along the c axis than 
in the ab plane, in contrast to our model prediction 
^ ^ab^ This suggests that other factors could con- 
tribute to the susceptibility anisotropy^i. For instance, 
anisotropic contributions are expected to arise from the 

Van Vleck susceptibility Xvv 7^ Xvv^ ^^'^ 
g- value anisotropy in the ground-state Kramers doublet, 
i-G-, gab 7^ 9ci due to lattice distortions lowering the oc- 
tahedral symmetry. These effects involve the inclusion of 
higher, excited spin-orbital, states. 

Nevertheless, we believe that the Kitaev-Heisenberg 
model describes some characteristic features of lay- 
ered iridates: for instance, recent susceptibility 
measurements^ on polycrystals of Li2lr03 and of 
Na2lr03 showed that the zero-temperature susceptibil- 
ity was about twice higher for the first (Li) compound 
than for the second. Assuming that both compounds can 
be described by the Kitaev-Heisenberg model and are in 
the stripe phase, the value of a for the Li compound 
should be significantly larger than for the Na compound 
- probably an > 0.6 if one compares the paramagnetic 
Curie temperatures of both compounds and follows the- 
oretical results of Ref. [s^ This is actually in agreement 
with the law x ^ jz^ which we find for any field di- 
rection (and thus for polycrystalline samples): since by 
assumption aNa > ccjv/s — 0.4, the observed ratio^ 
x{Li)/x{Na) ~ 2 would imply aLi > 0.7. This suggests 
that the Li-compound is quite close to the stripe-liquid 
transition. 

In the following, we shall confine our study to the 
Kitaev-Heisenberg model and explore the role of non- 
magnetic vacancies, unavoidably present in real com- 
pounds, and which, in general, can significantly alter the 
properties of frustrated magnetic systems^i^. 



III. COMBINED EFFECT OF A SINGLE 
VACANCY AND A MAGNETIC FIELD ON SPIN 
CORRELATIONS 

We now study the effect of a single non-magnetic va- 
cancy on the magnetic response of the system, especially 
at low-field. As previously, we will address separately the 
different phases of the model. 



A. Anomalous magnetization around a vacancy in 
the liquid phase 

We start by considering the liquid phase, and choose 
the field direction along the [001] axis, for which features 
characteristic of the Kitaev model are the clearest. First, 
we show that the combination of an infinitesimal field and 
a spin vacancy induces a peculiar magnetization pattern 
around the vacancy. The relevant quantity is the induced 




8a': h j= 0.05, a=l 50j': h ^= e, a= 0.9 



FIG. 6: (color online) (a, b) Correlation functions (afaj) 
with respect to site i indicated by a black circle, (c, d) In- 
duced magnetization patterns Saj due to the presence of a 
vacancy, on the N — 24 cluster. Circle areas are propor- 
tional to the quantities shown. Parameters are respectively 
{a,h^) = (1,0.05) in (a) and (c), while {a,h^) = (0.9,0) in 
(b) and (d) [a small symmetry-breaking field = e = 5.10~^ 
was added in (d)]. 



magnetization: 

which gives, for a given field hz and site j, the change of 
the magnetization at this site relative to its value (mo) in 
the vacancy-free system. We relate the vacancy-induced 
magnetization pattern to the correlation function (erf (7|) 
of the vacancy-free system. The most striking feature 
of the magnetization patterns shown in Fig. [Bfc.d) is the 
large magnetization on the z-neighbor of the vacancy (re- 
ferred to as site k in Fig. [7] from now on). The magneti- 
zation is significantly larger there than on any other site 
~ and only on this site does it remain finite at vanish- 
ing field in the Kitaev model: ((t|) Sj^k^nN for a = 1 
and hz ^ {m,N 0.48 for N = 24). Actually, in 
the thermodynamic limit (TL) the magnetization (both 
the total one and (tr^)) should behave with field asi^ 
hz\ii{hz)] our results are consistent with this prediction 
in the sense that m^v decreases with increasing iV, so 
that it should vanish in the TL as expected. As for the 
logarithmic field-dependence, we stress that this feature 
pertains to the TL, while for finite systems we can only 
notice weak non-linearities in the magnetization curves, 
discussed later in more detail. 

The magnetic response localized on the z-neighbor 
of the vacancy is clearly connected with the fact that 
the correlation function (cf crj) is non-zero only between 
nearest-neighbors sharing a z-bon d^^'^^ : if one spin is 
removed at site v, the term crfcr^ obviously disappears 
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from the Hamiltonian and thus the spin k is easily polar- 
izable by the z-oriented field (one can also interpret this 
behavior in terms of dangling Majorana fermion^i^) . 

Yet, we already saw in Section fll Al that a small field 
induces (erf (t|) correlations beyond nearest neighbors, it 
makes sense to focus in more detail on the magnetization 
pattern around a vacancy at small finite hz- Fig. EKc) 
shows the variation of magnetization (5cr| due to the pres- 
ence of the vacancy. Indeed, the spatial variations of this 
quantity presents similarities with those of correlations 
in the vacancy-free system: Sa^ is positive everywhere 
and decreases with increasing distance to the vacancy. 
Moreover, the spatial oscillations of correlation function 
appear in the magnetization pattern. Indeed, if one com- 
pares values of (5(t| at relative position either Hi + n2 or 
2ni — n2 from the vacancy [see values 0.09 and 0.013 re- 
spectively in Fig.inKc)], their respective ratio is compara- 
ble to that of corresponding correlations in the vacancy- 
free system. 

While (cr|cr|) correlations scale as hi, the magneti- 
zation enhancement JcrJ on further neighbors of the 
vacancy grows approximately linearly with h^; in con- 
trast, for J ~ k, it decreases from tun at /i = towards 
small negative values in the polarized phase (see 5a^ in 
Fig. [7]). Although the response observed here is obtained 
for a finite cluster size, we may expect that in a large 
system, the spin at site k acquires a large magnetization 
due to the combined effect of the neighboring vacancy 
and the small field, and in turn polarizes its other 
neighbors in a way that the (5cr| reproduces the above 
discussed correlation function pattern (Fig. [Sl[a)). 

If one considers the effect of a vacancy on the magnetic 
response of the system in a field hz, a relevant quantity to 
consider is 5M^{hz) = M^[N ^ 1, hz] - (iV - l)mo(/i^), 
with mo the magnetization per site of the vacancy-free 
system, and M^'[N — 1, hz] is the total magnetization of 
the system with — 1 spins plus a vacant site. 

This quantity, shown in Fig. [71 vanishes in the infi- 
nite field limit (where all spins are polarized), but also at 
small field in the TL. In finite systems, the previously dis- 
cussed finiteness of m^v goes along with a deviation from 
the expected logarithmic behavior of M^; yet, we notice 
on large enough systems (A^ = 24) a non-linearity of the 
magnetization curve in the small field range hz <SC hj^/p, 

the slope '^'^^J^ decreasing moderately with increasing 
small hz- But when hz becomes close to the value hj^/p 
of the liquid/polarized transition for the undoped sys- 
tem, the slope '^'^^^^ ^ increases strongly, so that (5M^ 
peaks near hz — 0.12(2). The exact position of the peak 
depends slightly on the cluster size, and the peak ampli- 
tude grows with increasing system size (roughly linearly: 
compare amplitudes for A^ = 16 and A^ = 24 on Fig. [7|). 
The existence of this peak can be understood in the fol- 
lowing way: not only does the presence of a vacancy cause 
an increase of the magnetization in the liquid phase, but 
it shifts slightly the liquid/polarized transition towards 




FIG. 7: Spin-vacancy induced magnetization SM'^ as a func- 
tion of magnetic field hz, for different a and cluster sizes. 
Shown also is a partial contribution (open circles) of the spin 
at site k, being the z-neighbor of the vacancy, for a = 1 and 
TV = 24. 



smaller fields h'j^^p < hj^/p (the difference between both 
fields being likely proportional to the vacancy concen- 
tration in a diluted limit), and the peak signals a field 
range h'j^^p < hz < hp/ p. This behavior can loosely 
remind of spin polarons formed around mobile holes in 
a paramagnetic background, so that the temperature of 
transition to the ferromagnetic (FM) phase shifts up with 
increasing hole density^. Here, instead of temperature 
the parameter driving the transition is the magnetic field, 
and moreover a vacancy is here not mobile, but its effect 
is as well to polarize its surroundings more and more 
when approaching the FM phase; consequently, at small 
finite vacancy concentration, the transition field should 
be slightly smaller than in the vacancy-free case. 

An analogy between the spin correlations in the 
vacancy-free system and the magnetization pattern 
around a vacancy can also be noticed in the case where 
the Kitaev Hamiltonian is perturbed by small Heisenberg 
interactions [see Fig. IHJb) and Fig. HKd)]. Note that in 
Fig. HUd) , we also included an infinitesimal field hz , the 
effect of which is only to break the Kramers degeneracy 
related to the odd number of spins. There, both the cor- 
relations and the vacancy-enhanced magnetizations (5cr| 
decrease rapidly with distance to the vacancy. Moreover 
the sign of these quantities alternates consistently with 
the pattern of z-oriented stripes shown in inset of Fig. HI 
and their decrease is much faster when going perpendic- 
ularly to stripes than along stripes: for the magnetiza- 
tions (cr|) (or equivalently 5a j, since mo — 0), one can 
again compare in Fig. [HJd) the absolute values at po- 
sitions 2ni - n2 (0.004) and ni + n2 (0.014) from the 
vacancy. 

The staggered z-bond magnetization, which measures 
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the z-stripe order parameter, 

^'(^^) = ^ E ^'"^-''(^hA + <^hB). (10) 

r 

is shown in Fig. |8l It also evidences that, in the hquid 
phase, the short-range stripe pattern is induced around 
the vacancy by smaU Heisenberg interactions. Starting 
from the a ~ 1 hmit, the magnetizations at further neigh- 
bors increase (in absolute values) significantly with in- 
creasing amplitude of Heisenberg interactions, i.e. with 
decreasing a. We expect that the extensive quantity 
Nm^{Qz) diverges in the TL for a — )■ as/L- Moreover, 
we note that the vacancy-induced magnetization SM^, 
for fixed small field h^, is reduced by the presence of 
Heisenberg interactions. It is significantly smaller for 
a = 0.9 than for a = 1, as the tendency to stripe or- 
der around the vacancy, at small field, competes with 
the aligning effect of the field. Consequently, the peak 
in 5A'F-{hz) shifts slightly towards higher fields with de- 
creasing a: the effect of Heisenberg interactions is to 
increase the position h'j^jp of transition to the polarized 

phase (as for /il/f in the vacancy- free case, see Fig. 

We briefly comment here on the fact that in the po- 
larized phase (and unlike in small field), the vacancy- 
induced magnetization change (5a^ at the z-neighbor is 
negative. This can be understood with the perturbative 
treatment of transverse interactions described in Section 
III At but here, to evaluate the quantum correction to the 
magnetization at a given site, we consider the off-diagonal 
(perturbing) couplings only on neighboring bonds - see 

Ref. [3^ for a similar approach. Within this approach, 

2 

the magnetization at this site is (cr^) ~ 1 — +2a~iy ' 
and the predicted value for 5a^ is of the same order as 
seen in Fig. [7] and it decreases in a similar way as in- 
creases, sufficiently far away from the transition field to 
either the liquid or the stripe phase. Qualitatively, one 
can interpret this by the fact that at high field, fluctua- 
tions from a fully-polarized state are effectively stronger 
at the vacancy's z-neighbor than elsewhere, since the ab- 
sence of ferromagnetic interaction on the missing z-bond 
makes it easier to flip this spin than other spins. With a 
similar reasoning one can explain that, in the same field 
range, the magnetization at other nearest neighbors of 
the impurity is enhanced compared to toq- 

If, instead of a field along an easy axis, a field in the 
[111] direction is applied, one can assumei^ that the mag- 
netic response is approximately a linear combination of 
responses to each field component: thus each neighbor of 
a vacancy acquires an anomalously large magnetization 
- the neighbor spin at site j, that would be connected 
to the missing spin by a 7-bond, has a magnetization 
(crj) ~ h^\n{h~f), and in turn brings a smaller polariza- 
tion, also parallel to h^, to other spins in its vicinity - 
mainly on neighboring sites which are close to the spatial 
axis in the prolongation of the missing 7-bond. Hence, 
three domains should coexist in the system, each with 
a different magnetization direction. Note that a vector 



'^stripe — {Qx),'my{Qy),'m^{Qz)] can be defined and 
used to characterize the polarization of spins around the 
vacancy. Qualitatively, rhstripe should point in a direc- 
tion closer to the easy axis 7 for which is the largest 
than to other easy axes; but it can have three non-zero 
components simultaneously. 



B. Stripe orientation in a small field and 
anisotropy of magnetic response 

In the stripe-ordered phase, among the 6 allowed stripe 
patterns (3 possible stripe directions times 2 spin ori- 
entations), a vacancy in combination with an infinites- 
imal field hz > Q selects the z-stripe pattern for which 
= 1 up to small quantum corrections. This selec- 
tion is due to the necessary imbalance between the num- 
ber of spins respectively aligned and antialigned with 
the field (see inset of Fig. This state is stable ver- 
sus increasing ft,^, until the field reaches a critical value 
(or anisotropy field han), where canted stripe patterns 
- those with stripe directions along x or y - become 
more favorable. This effect is similar to the local ro- 
tation of the staggered magnetization close to a vacancy 
in SU{2) antiferromagnets2^. Here, the main difference 
is that there is only a finite number of stripe patterns 
which cannot be continuously connected to each other 
without going out of the ground state manifold; conse- 
quently, the system has to choose between allowed stripe 
patterns, tilted or not by the field, and rhstripe is forced 
to point exactly along one easy axis 7, corresponding to 
(one of) the most favorable pattcrn(s). To understand 
the existence of a finite anisotropy field, one has to com- 
pare the energy of the z-stripe pattern {untiltable^ in a 
[001]-field) to that of the tiltable x- and y-stripe patterns. 
The former, thanks to its total magnetization ~ -1-1, 
acquires locally an energy gain ~ — /i^ with respect to the 
zero-field value but it has zero susceptibility [x^ = 0) 
thus no contribution cx NKj^ to its energy; in contrary, 
the field lowers the energy of the latter patterns by an 

extensive quantity 5E ~ —{N — 1) g(-]^^^^ ■ This compar- 
ison gives han = — 0-10 f^i' = 0.7, not far from 
the position of the jump in M'^ihz) in Fig. [H 

In the vicinity of the stripe-liquid transition, the mag- 
netization pattern induced by a small [001]-field {\hz\ <C 
0.1) around a vacancy can be of special interest: assum- 
ing that vacancies are diluted enough such that their mu- 
tual interactions can be neglected, the patterns expected 
on both sides of the transition are clearly distinct from 
each other. In the stripe phase, all neighboring spins of 
the vacancy have a large, almost field-independent mag- 
netization, with two of them aligned with the field and 
the third one antialigned with h. In contrast, in the liq- 
uid phase only one nearest neighbor spin is significantly 
polarized (with a strongly field-dependent polarization) 
and aligned with the field; and the two others arc weakly 
polarized opposite to the field. Thus, if one would have 
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FIG. 8: (color online) The 2-stripe order parameter m^(Qz) 
computed on the N = 24 cluster, as a function of a at different 
small fields hz. Right-inset: The magnetization pattern (a^) 
around the vacancy for a = 0.7 and = 10~^. Left-inset; 
Total magnetization as a function of a field applied in 
the liquid (a = 0.9) and in the stripe (a = 0.7) phases, with 
a magnetization jump at han — 0.062. 



a local probe to determine the (site-resolved) magnetiza- 
tion in the neighborhood of a vacancy - this probe could 
be '^^Na NMR for instance - the measured pattern could 
indicate on which side of the transition the system is. 

We have also seen that a vacancy can block the spin 
canting mechanism in the stripe phase, under a small 
field in (or close to) the [001] axis. This phenomenon is 
absent for a field in the [111] direction, where all stripe 
patterns remain equivalent and continuously tiltable by 
the field. Consequently, in a finite system, the situ- 
ation concerning the susceptibility anisotropy at zero- 
field: = < x'^^^': is now opposite to the previ- 
ously discussed case without vacancies where we found 
-^[111] < Concerning now the thermodynamic limit 
with a finite concentration n„ of vacancies, if one assumes 
that Hy is large enough to infiuence the magnetic response 
at low temperatures but small enough so that vacancies 
behave independently from each other, a similar effect 
could account for the experimentally observed x"'^ < x'^- 
However, given that in a small field, within the Kitaev 
model, the magnetization enhancement caused by two 
vacancies residing on the same sublattice is much larger 
than twice the magnetization enhancement caused by a 
single vacancy^, one should take care about cooperative 
effects between the vacancies. 



IV. INTERACTION BETWEEN VACANCIES: 
BROKEN ORIENTATIONAL SYMMETRY AND 
MAGNETIZATION PROPERTIES 

We turn now to the case where several vacancies are 
present in the system; we first consider their effect on 



the magnetization curves; then, we will focus on their 
effective interaction in the stripe phase and show how a 
vacancy pair can select an ordered pattern. 



A. Comparison of magnetic response with vacancy 
pairs in various phases 

It is instructive to compare the effect of vacancies or 
pairs of vacancies on magnetization curves, for given field 
direction and system size. To do this comparison we 
choose a [001]-field and consider the periodic N = 24 
cluster, but discuss as well implications for the ther- 
modynamic limit. Magnetization curves M'{hz) with- 
out or with vacancies, in different phases of the Kitaev- 
Heisenberg model, are shown in Fig. [9] 

In the Neel phase, although the total magnetization 
= J2i{'^i) is not a conserved quantity in presence of 
small Kitaev interactions (a <^ ajy/s), the magnetiza- 
tion curve M^{hz) on finite clusters displays a succession 
of plateaux (almost fiat for a = 0.1), as seen in Fig.[9ja). 
They are separated by jumps A(M^) ~ +2 correspond- 
ing to a fiip of one spin 1/2. At fixed a the step-like 
increase of the magnetization per site uiQ^hz) disappears 
only in the TL, the steps becoming smaller and closer to 
each other as A'^ increases. Moreover the steps disappear, 
at fixed cluster size, as a gets closer to a^/s- 

Within this phase, one is tempted to define an effec- 
tive spin 1/2, Seff, associated to the vacancy, such that 
M^{hz — > 0) = +1 corresponds to Seff pointing along 
+z; with two vacancies, their effective •S'^j^'s add up if va- 
cancies are on the same sublattice (e.g., n.n.n. case) and 
cancel each other if on opposite sublattices (n.n. case). 
From the behavior of M^ih^ — >■ 0) one obtains similar 
conclusions concerning effective couplings between other 
components of Seff- One could therefore think that an 
effective interaction between two vacancy spins, if it can 
be defined, should be 5'C/(2)-invariant, and possibly of 
Heisenberg type Slfj-S^j^; yet, the gapless nature of the 
undoped system and the continuous behaviour of mag- 
netization curves in the TL make it unclear how to for- 
mulate the effect of vacancies in terms of effective spins. 

In the liquid phase, in contrast, the magnetization does 
not show any plateau, but vacancies cause peculiar fea- 
tures at small fields: We have already seen that, com- 
pared to the zero- vacancy curve, the slope of M^{hz) at 
small field is significantly increased by the presence of 
one vacancy in the system. If two vacancies are present, 
their effect depends strongly on their respective sublat- 
tices (see Fig. E^c)): if they are on opposite sublattices 
the M^{hz) curve is, as in the vacancy-free case, linear 
for small hz [and x^i^z — >■ 0) is only slightly modified]; 
but if they are on the same sublattice, this slope is greatly 
increased - for a large range of the vacancy-induced 
increase in is more than twice that of the single- 
impurity case. Again, even though we cannot identify a 
logarithmic behavior x^i^z) ^ l/(^z In'^^^ h^) as derived 
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FIG. 9: (color online) Magnetization curves M'^{hz) on the 
cluster N — 24 with 0, 1, or 2 vacancies. The labels 
PI, ...,P4 correspond to relative positions of 2 vacancies as 
shown in panel (c). Panels (a),(b), and (c) correspond to 
a — 0.1, 0.7, 0.9, respectively. 




FIG. 10: (color online) Structure factors (a) S'^{Qz) and (b) 
S^(Qx) a± h = Q calculated for various respective positions of 
2 vacancies PI, P2, P3, indicated in the insets. Also shown is 
the structure factor of the vacancy- free system (open circles), 
for which S^(Qz) = S^{Qx). All results are obtained on the 
Af = 24 cluster. 



B. Selection of stripe orientation by vacancy pairs 



in Ref. [T^ for a = 1, we see that vacancies, as soon as 
they are not equally distributed on both sublattices, act 
hke partially-unbound moments causing non-hnearities in 
the magnetization properties at low field. 

In the stripe phase, as in the liquid phase, Kitaev in- 
teractions are strong enough to prohibit the existence 
of plateaux in the magnetization curves M^(hz), except 
at small fields in presence of vacancies: with one va- 
cancy, the previously discussed = 1 plateau is again 
suggestive of an effective spin 1/2 of the vacancy; and 
with two vacancies on the same sublattice, the Af^ ~ 2 
plateau seen on Fig.[9l^b) (pair PI) extends up to a field 
han2 < 0.22(1), that is more than twice the anisotropy 
field han for a single vacancy. This indicates a coopera- 
tive effect of both vacancies, which stabilize the z-stripe 
pattern more than if they behaved independently from 
each other. To describe this cooperative behavior, we 
will, similarly as in Section IIII Bl analyze the effect of 
vacancies on the energies of the various stripe patterns. 



We focus now on the stripe-ordered phase, and will 
explore how at h — 0, a pair of defects can, depending 
on their relative position, favor or disfavor one of the 
three stripe directions. In terms of symmetries, this can 
be understood as a breaking of Z3 symmetry related to 
the spatial orientation of stripes, while the Z2 symmetry 
((7i — >■ — (7i) remains unbroken. Remarkably, as shown in 
Fig. [ini the selected stripe pattern for a given vacancy 
pair depends on the coupling constant a. 
This effect is easy to understand for nearest-neighbor de- 
fects on a z— bond (configuration PI: see Fig fTU)) : classi- 
cally, the 2— stripe pattern has an energy lower by 

5E^ = 2 - 4a (11) 

than the one of other stripe patterns. Indeed, in the 
z-pattern only one FM bond - here a z-bond, with a 
classical contribution to the energy 1 — 3a - is lost due 
to vacancies, instead of two for other patterns (e.g. two 
x-bonds, in a a;-stripe pattern), and for a > 0.5 FM in- 
teractions contribute more than AF ones (whose classical 
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contribution to the energy is a — 1). For a < 0.5 the sit- 
uation is reversed, which explains that S^{Qz) becomes 
small and S'^iQ^) = ^^(Qy) ~ 1/2. 

For second neighbors (case P2 on Fig. [TU]) , stripe pat- 
terns are classically degenerate, but this degeneracy is 
lifted by quantum fluctuations. We evaluate in second 
order perturbation theory the contribution to the energy 
of various bonds in which at least one site is a neighbor 
of the vacancy; actually only AF bonds give a contribu- 
tion in this approximation. For instance fluctuations on 
a bond belonging to the same hexagon as the two impuri- 
ties yield a contribution -4(l-2Q;)2/[4(l-a) + 2(3a-l)] 
to the energy of the z-stripe state. Collecting contribu- 
tions from all relevant bonds, the energy of the z-stripe 
state (relative to that of other stripe orientations) is: 



6E^ = 2(1 - 2ay 



'5a — 1 2a 1 + a 



(12) 



which is positive for a in the vicinity of 1/2, and changes 
sign at a2 — 0.651. This explains the drop of the z-stripe 
structure factor (dotted curve on Fig. [TUTa)). from large 
values for a £ [0.64(2); 0.8], to almost zero for smaller a. 
The smaller values of S^{Qx) — 0.3 in the vicinity of a = 
0.5 are related to the smaller energy scale separating the 
different stripe patterns: the above expression, with its 
prefactor (1 — 2a)^, indicates that the 3 patterns become 
degenerate at a = 0.5. 

A similar energy comparison between different pat- 
terns for third neighbor vacancy pairs (configuration P3: 
see Fig. [TO)) predicts a change of stripe orientation at 
as = 0.6 between z— stripes (as < a < uls) and x- 
stripes (asN < a < as), following from the perturbative 
result: 



5EI 



(1 - 2a)2(3- 5a) 



(13) 



whereas this orientational change is seen on = 24 at a 
slightly smaller value a^ ~ 0.53(1). Notice that here (and 
more generally for pairs of vacancies at further neighbor 
sites which are not related to each other by any reflec- 
tion symmetry of the honeycomb lattice) the three stripe 
orientations are all non-equivalent to each other. 

Finally we note that at zero-field, the conclusions 
above hold also for pairs of neighbor vacancies with other 
relative positions: for instance with next-nearest neigh- 
bor vacancies separated by a x- and a z-bond and a > a2 , 
according to the previous reasoning the favored pattern 
is a j/-stripe pattern. 

In presence of a small magnetic field, the Z2 symmetry 
characterizing each stripe orientation (e.g., z-stripes) is 
broken. Moreover, similar effects occur as in the previ- 
ously discussed single vacancy case: if the most favor- 
able pattern in an infinitesimal [001] field (e.g., z-stripe 
pattern around a vacancy pair in configuration PI, for 
a > 0.5) is untiltable for this field direction, as the field 
amplitude gets larger than a critical value han2 this 
pattern becomes less favorable than other patterns which 



are tillable. In the limit of large systems A'' ^ 1, a clas- 
sical estimate of han2 gives a result 4(2a — 1)(1 — a)/\/N 
which, for a = 0.7 and A^ = 24, has a value han2 — 0.28 
comparable to the observed position of the magnetization 
jump. 

The relevant quantity for the discussion of either low 
temperature properties or magnetic ordering at small va- 
cancy concentrations is the energy scale 5E^^, namely the 
energy which favors a particular stripe pattern in the 
vicinity of a vacancy pair - n standing for the number 
of bonds separating both vacancies. Far away from the 
stripe/liquid transition (say, for a < 0.7) 5E^ decreases 
rapidly with n (e.g., \5E^\ < \8El\/2 for n>2). Closer 
to as/Lj larger quantum fluctuations make the stripe- 
orientation selection by 3'"'' neighbor vacancy pairs as 
efficient as by nearest neighbor vacancy pairs. 

Several consequences can result from this: (i) In the 
simpler case with only two vacancies, the energy lower- 
ing of some stripe patterns for a ~ Q-s/l also stabilizes 
these patterns with respect to the spin liquid, which leads 
to a slight shift of the stripe/liquid transition towards 
larger a - possibly 0.82 to 0.85 depending on the relative 
position of vacancies, as can be seen in Fig. [TUTa). (ii) 
More importantly, if one focuses on systems with sub- 
stantial vacancy concentrations, say, ^ 10%, neigh- 
boring vacancy pairs as described in Fig. [10] are abun- 
dant enough to select locally stripe orientations that will 
frustrate each other at the global scale, (iii) As a result, 
the system could display a glassy behavior, especially in 
the vicinity of the stripe-liquid transition. The resulting 
glassy phase would consist of stripe-ordered nanoscale 
domains, with a vacancy pair in each domain, selecting 
its stripe orientation. In order to perform a modeling 
of such a system and to reproduce a possible spin-glass 
behavior, one should obviously need to consider systems 
of (at least) hundreds of Ir atoms, and to take into ac- 
count the energy cost of domain walls between regions of 
different stripe orientations. Moreover, obviously not all 
vacancies are grouped into distinct first, second or third 
nearest neighbor vacancy pairs; in particular the effect of 
unpaired vacancies, even at zero external field, could be 
subtle, although one can imagine single vacancies pinning 
walls between domains of different stripe orientation. 



V. SUMMARY AND CONCLUDING 
REMARKS 

We have considered here several important aspects of 
the Kitaev-Heisenbcrg model, which, despite its simple 
structure, displays in its T = phase diagram three dis- 
tinct phases, two of which are antiferromagnetically or- 
dered while the third one is a spin liquid. We employed 
numerical methods (exact diagonalization) for our inves- 
tigation of the properties of the model in a magnetic field 
and the effect of vacancies on the magnetic properties. 
The numerical methods were complemented by classical 
and perturbative analysis particularly in the frame of the 
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stripe phase. 

We showed how the ground state phases evolve in a 
magnetic field oriented along the [001] direction (one of 
the local easy axes of the Kitaev-Heisenberg Hamilto- 
nian). We found that the three phases are robust to 
small fields. Next, we analyzed the influence of the field 
direction. The magnetic response is qualitatively distinct 
in the different phases. In the stripe phase, the magne- 
tization can be described by a spin canting mechanism, 
which also provides an understanding of the anisotropy 
in the numerically obtained susceptibilities. In the liquid 
phase, we compared the effect of small magnetic fields 
and of small Heisenberg interactions, each considered as 
a perturbation of the Kitaev Hamiltonian. Such per- 
turbations are shown to induce specific patterns in the 
spin-spin correlations beyond nearest neighbors, while at 
larger fields a transition into the spin-polarized phase oc- 
curs. We determined the phase diagram of the model in 
a [001]-field, with good agreement between numerical re- 
sults and analytical estimates. For the latter analysis, we 
took into account both the spin canting and the quan- 
tum fluctuations in the ordered and polarized phases. 
This phase diagram has similarities with its counterpart 
for a [111] field orientation^^, however, a major difference 
between them is that for a [001]-field we do not see any 
quantum critical point in the vicinity of the stripe-liquid 
phase transition. We presume that this is related to the 
absence of a field- induced gap in the low-energy spectrum 
of the liquid phase for a [001]-fieldii. 

Subsequently, we focused on the effect of non-magnetic 
vacancies on the spin correlations and the magnetic re- 
sponse of the model. For a single vacancy in the liq- 
uid phase, we found that the spatially anisotropic mag- 
netization pattern around the vacancy is related to the 
anisotropy of spin correlations in the vacancy-free sys- 
tem. Up to finite-size effects inherent to our approach, 
these results are also in agreement with predictions of a 
non-linear response to a small [001]-field, whose effect is 
mostly seen on one specific neighbor of the vacancy. In 
the stripe phase, a vacancy coupled to a small [001]-field 
breaks both the Z2 (time reversal) and the Z3 (orienta- 
tional) symmetry characterizing the ground state mani- 
fold; a single stripe pattern is selected, which is not con- 
tinuously tillable by the field in contrast to the vacancy- 
free situation. 

Having in mind the description of systems at finite va- 
cancy concentration, we considered the interplay of two 
vacancies sitting either on nearest or further neighbor 
sites. Again the magnetic response to an easy-axis field 
depends strongly on the nature of the phase, but it also 
depends sensitively on the relative position of vacancies. 
In the liquid phase a non-linear behavior, much stronger 
than for the single-impurity case, is seen if the two va- 
cancies are on the same sublattice, while the response 
remains linear if they reside on different sublattices. In 
the stripe phase, vacancy pairs select a specific stripe ori- 
entation even at zero field. This selection mechanism, de- 
pending on the value of coupling parameter a, can be well 



understood by a pcrturbative approach evaluating locally 
the effect of quantum fluctuations. The response to a 
small field depends on whether the selected stripe pat- 
tern is tillable or not for this field direction2£. Although 
the strongest selection effect is caused by nearest neigh- 
bor vacancy pairs for a wide range of a, further neighbor 
pairs have to be considered as well in the vicinity of the 
stripe-liquid transition. Therefore we conclude that the 
intrinsic randomness of vacancies, in combination with 
the selection mechanism induced by vacancy pairs, may 
lead to spin-glass behavior in the stripe regime. There, a 
vacancy concentration of a few percents should be enough 
for the system to show a glassy behavior, with a tendency 
to the formation of stripe-ordered nanodomains around 
vacancy pairs. 

Finally, we comment on the possible relevance of this 
work to layered iridates, which motivated the derivation 
of the Kitaev-Heisenberg model. The experimental data 
(magnetic susceptibility, heat capacit y^^i^" ) on Na2lr03 
indicate antiferromagnetism at low temperature and a 
substantial magnetic anisotropy. Yet, the angular de- 
pendence of the measured susceptibility is qualitatively 
opposite to that found for a vacancy-free model. We 
then showed that non-magnetic vacancies in the dilute 
limit, i.e., where interactions between vacancies are ir- 
relevant, can reverse locally this tendency by blocking 
the spin-canting for some field directions. The recent x- 
ray magnetic scattering data^i ruled out the Neel state, 
implying that the Kitaev interaction in Na2lr03 may 
dominate over a simple Heisenberg coupling. In addi- 
tion to the stripe phase intrinsic to the Kitaev-Heisenberg 
model, the so-called "zig-zag" spin order is also consis- 
tent with the data^i. An element in favor of a descrip- 
tion of layered iridates A2lr03 (A=Li,Na) by the Kitaev- 
Heisenberg model is the comparison of the low temper- 
ature susceptibilities of these compounds^. Using our 
result X oc jzr^ for the stripe phase, we found a rela- 
tion between a values in these two compounds, which 
indicates that the Li compound, as the Na-based one, 
is most probably in the stripe phase but much closer to 
the transition towards the spin liquid. This conclusion 
is in agreement with previous work^^iS estimating the 
relative a values from a comparison of the paramagnetic 
Curie temperatures. The application of uniaxial pressure 
on a single crystal of Li2lr03 might even drive the system 
into the liquid phase. 

To describe the experimental situation in layered iri- 
dates on a quantitative level, other factors such as lattice 
distortions have to be considered as well. In particu- 
lar, the lattice distortions may bring about substantial 
anisotropy in the ground state g-factors and in the Van 
Vleck contribution to the magnetic susceptibility, origi- 
nating from transitions to higher lying spin-orbital quar- 
tet split by non-cubic (tetragonal and/or trigonal) crystal 
fields. Furthermore, these distortions may lead to spin in- 
teractions not included in the Kitaev-Heisenberg model. 
For intermediate strength of spin-orbit coupling like in 
iridates, these effects might be essential for a quantita- 
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tive description of magnetic properties. Moreover, our 
study suggests that a glassy behavior could occur in the 
vicinity of the spin liquid phase; in this context, con- 
sidering the present model at larger scales, with a finite 
density of vacancies, could help to explore the effect of 
temperature. 
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